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Anomalous diffusion in nonhomogeneous media: Power spectral density of signals 
generated by time-subordinated nonlinear Langevin equations 
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Subdiffusive behavior of one-dimensional stochastic systems can be described by time- 
subordinated Langevin equations. The corresponding probability density satisfies the time-fractional 
Fokker-Planck equations. In the homogeneous systems the power spectral density of the signals gen¬ 
erated by such Langevin equations has power-law dependency on the frequency with the exponent 
smaller than 1. In this paper we consider nonhomogeneous systems and show that in such systems 
the power spectral density can have power-law behavior with the exponent equal to or larger than 
1 in a wide range of intermediate frequencies. 
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I. INTRODUCTION 

A number of experimental observations show that more complex diffusion processes in which the mean-square 
displacement is not proportional to the time t take place in various systems. A broad family of processes described 
by certain deviations from the classical Brownian linear time dependence of the centered second moment is called 
anomalous diffusion. Anomalous diffusion in one dimension is characterized by the occurrence of a mean square 
displacement of the form 


((Axn = -2^t“, (1) 

r(l -b a) 

which deviates from the linear Brownian dependence on time [T| . Eq. m introduce the anomalous diffusion coefficient 
Ka- Such a deviation from classical diffusive behavior can be observed in many systems and leads to many 

interestiM plwsical properties Applications of anomalous diffusion have been found in physics, chemistry and 
biology [l|, la, |6| . In general, anomalous diffusion occurs in complex structures exhibiting the presence of long-range 
correlations or memory effects [l|. In the physics of complex systems, anomalous transport properties and their 
description have attracted considerable interest starting with the pioneering papers of Montroll and his collaborators 
0 . 

An important subclass of anomalous diffusion processes constitute subdiffusion processes, characterized by the 
sublinear dependence with the power-law exponent in the range 0 < a < 1. In this situation no finite mean jump time 
At exists . Subdiffusion processes have been reported in condensed matter systems 0 , ecology @ , and biology 0 . 
Continuous time random walks (CTRWs) with on-site waiting-time distributions falling slowly as and lacking 

the first moment predicts a subdiffusive behavior and is a powerful tool to describe systems which display subdiffusion 
[3, Q. Starting from the generalized master equation or from the CTRW the fractional Fokker-Planck equation can 
be rigorously derived 0 [T^. Fractional Fokker-Planck equation provides a useful approach for the description of 
transport dynamics in complex systems which are governed by anomalous diffusion Q and nonexponential relaxation 
patterns [ul . It has been used to model dynamics of protein systems and for reactions occurring in disordered media 
ElIMl- Description equivalent to a fractional Fokker-Planck equation consist of a Markovian dynamics governed 
by an ordinary Langevin equation but proceeding in an auxiliary, operational time instead of the physical time M- 
This Markovian process is subordinated to the process defining the physical time; the subordinator introduces memory 
effects [ 2 ^ . Other approaches for the theoretical description of the subdiffusion use the generalized Langevin equation 
[2ll - l^ , fractional Brownian motion , or the Langevin equation with multiplicative noise [ 2 ^ . 

The traditional CTRW provides a homogeneous description of the medium. More complex situation is the diffusion 
in nonhomogeneous media, for example diffusion on fractals and multifractals [ 2 ^. Nonhomogeneous systems exhibit 
not only subdiffusion related to traps, but also enhanced diffusion can occur: for example, transport of interacting 
particles in a weakly disordered media is superdiffusive due to the disorder and subdiffusive without the disorder [T7I| . 
Anomalous diffusion in heterogeneous fractal medium has been considered in Ref. [ 2 ^ where it was proposed that 
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in one dimension the mean square displacement has the form ((Ax)^) ^ instead of Eq. ([IJ. Hetero gen eous 

fractional Fokker-Planck equation on heterogeneous fractal structure media has been investigated in Refs. [29l - [^ . 
In nonhomogeneous media the properties of a trap can reflect the medium structure, therefore in the description 
of transport in such a medium the waiting time should explicitly depend on the position. This dependence can be 
introduced by using the position-dependent subdiffusion exponents Another way is to consider position- 

dependent time subordinator [^ . 

In the homogeneous systems the power spectral density (PSD) of the signals generated by time-subordinated 
Langevin equations has power-law dependency S{f) ~ on the frequency as / —?> 0. [s^. Since 0 < a < 1, the 
power-law exponent 1 — a is smaller than 1. The purpose of this paper is to consider the PSD in nonhomogeneous 
systems exhibiting anomalous diffusion. We demonstrate, that in such systems the PSD can have power-law behavior 
with the exponent equal to or larger than 1 in a wide range of intermediate frequencies. 

The paper is organized as follows: In Sec. HIl we introduce the time-fractional Fokker-Planck equation describing 
subdiffusion in nonhomogeneous media. The expression for the power spectral density of the fluctuations of the 
diffusing particle in such a medium is obtained in Sec. IIIII In Sec. IIVI we consider a particular case of the time- 
fractional Fokker-Planck equation involving the coefficients with power-law dependence on the position. Numerical 
methods of solution are discussed in Sec. IVl Section rvTl summarizes our findings. 


II. TIME-FRACTIONAL FOKKER-PLANCK EQUATION FOR NONHOMOGENEOUS MEDIA 

In this Section we derive the time-fractional Fokker-Planck equation describing diffusion of a particle in nonhomoge¬ 
neous media. Usually the description of the anomalous diffusion is given by the CTRW theory assuming heavy-tailed 
waiting-time distributions between successive jumps of the diffusing particle. Here we use the method of the derivation 
that is similar to that outlined in Refs. [13,1^. We start with the Markovian process described by the Ito stochastic 
differential equation (SDE) 


dx{T) = a{x{T))dT -\- h{x{T))dW{ t) . (2) 

Here W{t) is the standard Brownian motion (Wiener process). The drift coefficient a{x) and the diffusion coefficient 
h{x) explicitly depend on the particle position x. This dependence on the position reflects the nonhomogeneity of a 
medium. Following Ref. we interpret the time r in Eq. ([2]) as an internal, operational time. Equation @ we 
consider together with an additional equation that relates the operational time r to the physical time t. The difference 
between physical time t and the operational time t occurs due to trapping of the diffusing particle. For the trapping 
processes that have distribution of the trapping times with power law tails, the physical time t = T (r) is given by the 
the strictly increasing a-stable Levy motion defined by the Laplace transform 

(g-feT(r)) ^ ^ (3) 

Here the parameter a takes the values from the interval 0 < a < 1. Thus the physical time t obeys the SDE 

dt{T) = dL°‘( t) , (4) 

where dL°^{T) stands for the increments of the strictly increasing a-stable Levy motion For such physical time 

t the operational time r is related to the physical time t via the inverse a-stable subordinator [s^, |4^ 

S{t) = inf{r : T(r) > t} . (5) 

The processes x{t) and S{t) are assumed to be independent. Equations (H]) and (|1]) define the subordinated process 
y(t) obtained by a random change of time 


y{t) = x{S{t)). (6) 

The process y(t) describes the diffusion of a particle in a medium with traps. 

We will derive the equation for the probability density function (PDF) of y. For the derivation we use the method 
of Laplace transform. The PDF Px{x, t) of the stochastic variable a; as a function of the operational time t obeys the 
Fokker-Planck equation corresponding to the Ito SDE ([2]) 


d 

—Px{x,t) = Lpp{x)Px{x,t) , 


( 7 ) 
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where Lpp (x) is the time-independent Fokker-Planck operator [4l| 

The Laplace transform of Eq. 0 is 

kPx{x,k) - Px{x,0) = LFp{x)Px{x,k). (9) 

Since the processes x{t) and S{t) are independent, the PDF of the random process x{S{t)) is given by 

P{x,t) = j Px;{x,T)Ps{T,t)dT . (10) 

Here Ps{t, t) is the PDF of the inverse a-stable subordinator S{t). From Eq. (IlOp it follows that the Laplace transform 
P{x, k) of the PDF P{x, t) is related to the Laplace transform Ps(r, k) of the inverse subordinator S{t): 

P{x,k)= J Px:{x,T)Ps{T,k)dT . (11) 


The Laplace transform Ps(t, k) of the inverse subordinator S(t) we obtain as follows: from the definition of the inverse 
subordinator ([51) we have Pr(S'(t) < r) = Pr(r(r) > t), therefore 

PsiT,t) = -^J*PT{t',T)dt'. ( 12 ) 

Here PtP^t) is the PDF of the strictly increasing a-stable Levy motion T(r). The PDF PtP^t) fulfills the scaling 
relation 


1 


Prityx) = —Pt —, 1 , 


t 


(13) 


since the strictly increasing a-stable Levy motion is 1/a self-similar [d^. Combining Eqs. and we obtain 


Psixyt) = —PrityT). 
ar 

Consequently, the Laplace transform of Ps(r, f) is equal to 

Ps{T,k) = . 


Here we used Eq. m for the Laplace transform of Ppp, t). 
Using Eqs. m and m we get 


P{Xyk) = k^-^Px{xyk°‘). 


(14) 


(15) 


(16) 


Acting with the operator LFp(a;) on Eq. (fTCl) we have 

P{x, k) = k~^Px{Xy 0) -b k~^LFp(x)P{Xy k). (17) 

The inverse Laplace transform of this equation yields 

1 P 

P{Xyt) = Px{XyO) +-— dt' {t - t'P~^Lpp{x)P{x,t'). (18) 

r(a) Jo 

Introducing the fractional Riemann-Liouville operator [d^ 

oA““/W = - dt' y 0 < a < 1 (19) 

‘ r(a) 7 o (i-t') 

we can write Eq. as 


P{Xy t) = Px{Xy 0) -b oDt °‘Lpp{x)P{Xy t) 


( 20 ) 
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By differentiating this equation with respect to time we get the time-fractional Fokker-Planck equation 

Ap,.,,) = . (21) 

where 

0<»<i (22) 

The operator qDI~°‘ is expressed via the convolution with a slowly decaying kernel, which is typical for memory effects 
in complex systems Equation dan is the equation describing the subdiffusion of particles in an inhomogeneous 
medium. This equation generalizes the previously obtained time-fractional Fokker-Planck equation with the position- 
independent diffusion coefficient. 


A. Position-dependent trapping time 

The properties of a trap in a nonhomogeneous medium can reflect the structure of the medium. In the description 
of the transport in such a medium the waiting time should explicitly depend on the position [s^ . Instead of Eq. 0 
we assume that the physical time t is related to the operational time t via the SDE 

dt(.r) = g{x{T))dL°‘{T). ( 23 ) 

Here the positive function g{x) is the intensity of random time and models the position of structures responsible for 
either trapping or accelerating the particle. Large values of g{x) corresponds to trapping of the particle, whereas 
small g{x) leads to the acceleration of diffusion. Similar equation has been used in Ref. [^. We interpret Eq. (1^ 
according to the ltd stochastic calculus: the values of x and t at operational time r are determined by events prior 
to the application of the stochastic force which acts only from time r to t -|- dr. This assumption leads to the 
decoupling of the changes of x and the changes of t occuring during an infinitesimal increment of the operational time 
dr. Note, that the increments of the strictly increasing a-stable Levy motion L'^{t) are characterized by long tails 
and thus only moments of order smaller than a are finite. 

For fixed particle postion x the coefficient g{x) in Eq. (1^ is constant and Eq. (1^ corresponds to the fractional 
Eokker-Planck equation 


t\x) = -oD^g{x)°‘P{t; t\x) . (24) 

This equation can be obtained by noting that from the definiton of the stricly increasing a-stable Levy motion ([3]) the 
Laplace transform of the PDF P{t;T\x) is P{k-,T\x) = exp{—T[g(x)fc]“}. Diferentiating this expression with respect 
to T and taking the inverse Laplace transform one rets Eq. (l24l) . Alternatively, one can obtain the fractional Fokker- 
Planck equation using the methods of Refs. [2^, |4^, The fractional derivative in the Fokker-Planck equation 
appears as a consequence of the increments of Levy a-stable motion in Eq. (1231) . 

Equations 0 and (133)) together define the subordinated process. However, now the processes x{t) and ^(t) are 
not independent and the derivation of the Fokker-Planck equation presented in previous subsection is not applicable. 
Nevertheless, we can show that also with position dependent trapping time the resulting equation has the form of 
Eq. (1311) . To do this let us consider the joint PDF Px^t{x,t;T) of the stochastic variables x and t. 

SDEs ([3]) and (1331) correspond to the following two-dimensional fractional Fokker-Planck equation: 

= LFp{x)Px^t - oDfg{x)'^Px^t ■ (25) 

or 

This equation is a combination of Eqs. © and (1241) . Two-dimensional fractional Fokker-Planck equation (1351) for the 
PDF of two stochastic variables x and t can be rigorously derived from the SDEs © and (1331) driven by Levy stable 
noises as in Refs. [1^ (the Gaussian noise in Eq. (I3|) is a particular case of a Levy stable noise with index of 

stability a = 2). 

The zero of the physical time t coincides with the zero of the operational time r, therefore, the initial condition 
for Eq. (1331) is Px,t{x,t]Q) = Px{x,0)S{t). In addition, since t is strictly increasing, we have a boundary condition 
Px,t{x,0-,T) = 0 when r > 0. The fractional Riemann-Liouville operator in Eq. (1351) we can write as = 

d — l 

■ 
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Now let us consider x and r as stochastic variables instead of x and t. Since the stochastic variable t is related to 
the operational time r via Eq. the joint PDF P^ ,^{x,T]t) of the stochastic variables x and r is related to the 
PDF according to the equation 


^x,T (^5 : i') 


= oDt ^9{x) 


'‘Px,t{x,t;T) . 


(26) 


This equation can be obtained by noting that the last term in Eq. (1251) contains derivative ^ and thus should be 
equal to —^Px,T- From Eq. (E51) if follows that 


Pxx. — oD 


1 —Q! 


5(a 


-Px 


(27) 


Using Eqs. (1^ and (07)) we obtain 


—Px,rix,T;t)=QDl °‘Lfp{x) 


gix)' 


-Px,r g^oD 


gix) 


-Px 


(28) 


The PDF Px^T has the initial condition Px^rixjT^O) = Px{x,0)S{t) and the boundary condition Px^rix,0;t) — 0. The 
PDF of the subordinated random process x{t) is P{x, t) = J Px,t{x, t ; t) dr . Integrating both sides of Eq. (l28l) we get 


where the new Fokker-Planck operator is 


-Pix,t) = oDl-^L'Mx)P, 






Here the new drift and the diffusion coefficient are 

a{x) 


a [xj = 


5 (a 


b'{x) = 


bjx) 

gix)'- 


(29) 


(30) 


(31) 


Thus position-dependent trapping leads to position-dependent coefficients in the time-fractional Fokker-Planck equa¬ 
tion, even if the initial SDE ([7|) has constant coefficients. Eq. (I29|) is the same as Eq. (EH) when g{x) is constant and 
does not depend on position. 


III. POWER SPECTRAL DENSITY AND TIME-FRACTIONAL FOKKER-PLANCK EQUATION 

In this Section we derive a general expression for the PSD of the fluctuations of the diffusing particle in nonhomo- 
geneous medium. The evolution of the PDF of particle position x is described by the time-fractional Fokker-Planck 
equation (EH)- For calculation of the spectrum we use the eigenfunction expansion of the Fokker-Planck operator Lpp- 
Method of eigenfunctions for solving of time-dependent fractional Fokker-Planck equation has been used in Ref. [d^. 
Spectrum of fluctuations when the diffusion coefficient is constant has been obtained in Ref. (3^ . Similar derivation 
of the spectrum for nonlinear SDE has been performed in (d^ . 

The eigenfunctions of the Fokker-Planck operator Lfp{x) are the solutions of the equation 

LFpix)Pxix) = -XPxix). (32) 

Here P\ix) are the eigenfunctions and A ^ 0 are the corresponding eigenvalues. The eigenfunctions obey the or¬ 
thonormality relation |d9l | 


J e^^-^Pxix)Pyix)dx = , 


(33) 


where 


$(x) = — lnPo(a;) 


(34) 


is the potential associated with the operator LFp(a:)- Here Poix) is the steady-state solution of Eq. (ITTl) . 
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We can write the time-dependent solution of the fractional Fokker-Planck equation (EU corresponding to a single 
eigenfunction as 


P{x,t) = Px{x)fx{t). 

Inserting into Eq. (ED we get that the function f{t) obeys the equation 

^/aW = -AoA'-“/aW 

with the initial condition /(O) = 1. The Laplace transform of this equation yields 

kMk) = l-Xk^-^h{k). 

The solution of Eq. ED is 

^ k + Xk^-^^ ■ 

The inverse Laplace transform is given in terms of the monotonically decreasing Mittag-Lefller function (47 

hit) = E^i-xr). 

The Mittag-Leffler function has a series expansion 


(35) 

(36) 

(37) 

(38) 

(39) 


Ea{z) = Ea^i{z) = ^ 


r(an 1) ■ 

n —0 ^ ' 


(40) 


The autocorrelation function can be calculated from the transition probability P{x^t\xo,0) (the conditional prob¬ 
ability that at time t the stochastic variable has value x with the condition that at time t = 0 it had the value 
xo): 


C{t) = J dx J dxo xoxPo{xo)P{x,t\xo:0) — 


dx xPq (x) 


(41) 


The transition probability is the solution of the Eokker-Planck equation (ICT) with the initial condition P(x, Ojxo, 0) = 
S(x — xq)- Expansion of the transition probability density in a series of the eigenfunctions has the form 


P(x,tlxo,0) = J^Px(x)e^(^°^Px(xo)E^(-Xn, (42) 

A 

where we used Eqs. ed and Inserting Eq. (l42ll into Eq. ED we get the expression for the autocorrelation 

function 


c(t) = j^x^xPo.(-xn. 

A>0 


(43) 


Here 




xP\(x) dx 


(44) 


is the first moment of the stochastic variable x evaluated with the A-th eigenfunction Px(x). Such an expression for 
the autocorrelation function has been obtained in Ref. [s^ ■ 

According to Wiener-Khintchine relations, the power spectral density is related to the autocorrelation function: 

^OO 

S{f) = A / C (t) cos{ujt) dt, (45) 

Jo 

where w = 27r/. Using Eq. (I43II we obtain 


S{f)=AY,Xl 

\>o 



Ea{ — Xt°‘) cos(wt) dt 


(46) 
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The integral can be calculated by noting that the Laplace transform of Ea{—Xt‘^) is given by Eq. (l38l) . We obtain the 
desired expression for the PSD 


^(/) = 4- 


sin (fa) 


,1-Q! 


^ + 2Aw“ cos (fa) 




(47) 


Eq. (|47l) becomes the usual expression for the PSD when a —>■ 1. Similar expression for the spectrum has been 
obtained in Ref. [s^ . 

Eor small frequencies lu Aj^' we can neglect the frequency when it appears together with the eigenvalues A. Here 
Ai is the smallest eigenvalue larger than zero. Thus for small frequencies Eq. (SZl) approximately is 


S{f) 


sm(fa) ^ 

,1-a A 




X 


(48) 


We obtain that for small frequencies the PSD has a power-law dependency on the frequency S{f) ^ However, 

the power-law exponent is always smaller than 1, since 0 < a < 1. It is not possible to get pure 1// spectrum this 
way. In the next Section we show that it is possible to get larger power-law exponents in the PSD in a wide range of 
intermediate frequencies when the diffusion coefficient is not constant and depends on x. 


IV. TIME-FRACTIONAL FOKKER-PLANCK EQUATION WITH POWER-LAW COEFFICIENTS 


In this Section we consider a particular case of the time-fractional Fokker-Planck equation We assume that 

the diffusion coefficient has a power-law dependence on the particle position x and Eq. (EB takes the form 

|p(x,t) = {(^ -v)§-^ [x^^P{x,t)\ I . (49) 

Here rj is the power-law exponent of the multiplicative noise in Eq. (ED and V defines the behavior of the steady-state 
PDF Po{x). Eq. (HHl) should be considered together with the boundary conditions that restrict the stochastic variable 
X to the positive values. 

The steady-state PDF Po{x) obtained from Eq. ((4^ has a power-law form 

Po{x) ^ X-'. (50) 


For ly' > I the PDF Po{x) diverges as a; —>■ 0, thus the diffusion should be restricted at least from the side of small 
values. This can be done by introducing an additional potential that becomes large only when x acquires values 
outside of the interval [a;i„in, a^max] into the drift term of Eq. (|4^ . The simplest choice is the reflective boundaries at 

^ ~ a^min nnd X = X^nsiyi- 


fast electrons in a hot plasma 
nonlinear SDEs proposed in Refs. 


The power-law form of the diffusion coefficient is natural for systems exhibiting self-similarity, for example disordered 
materials, and has been used to describe diffusion on fractals [10, HH, turbulent two-particle diffusion, transport of 

10, ^ I ■ Equation (|4^ is a generalization of the Fokker-Planck equation resulting form 
54, HO] • Such nonlinear SDEs generate signals haviM 1/ f spectrum in a wide range 
of frequencies and have been used to describe signals in socio-economical systems |56l . l57l | and Brownian motion in 
inhomogeneous media [10]. 

In Ref. [iOj an approximate expression for the first moment Xx has been obtained for the Fokker-Planck operator 
appearing in Eq. ((4^ assuming reflective boundaries at Xmin = 1 and Xmax = C, C ^ 1- According to the results of 
Ref. S 


Aa 


CA 1 


(51) 


where 


jl-Tyj v-1 
- 




p = 


1^’ 


v — 3 


c\ = 


/3i = 1 + 


2 (^- 1 ) ■ 


(52) 
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The parameters Zmin and Zmax depend on the boundaries a^min and Xmax- When pzmax ^ 1, replacing summation by 
integration in Eq. (H71) we obtain the expression for the PSD 


S{f) « 4 


sin 





l — a 


A 

A2+a;2“ + 2Au;“ cos (fa) 


XlD{X) dX 


The density of eigenvalues D{X) has been estimated as 


El 


D{X) 


1 

Ta' 


Using Eqs. m and (l54)) we get 


S{f) ^ 4 


sin (f g) 

^l+a(/3i-l) 



1 

U^l-1 


1 

+ 1 + 2m cos (fa)) 


du 


(53) 


(54) 


(55) 


Here the upper range of integration is limited because X\ becomes small when pzmin » 1 0 . When < a;“ <C 
•^inhi ^-nd 0 < /3i < 2 then we can approximate the lower limit of integration by 0 and the upper limit by oo. In this 
case the PSD depends on the frequency as S{f) ~ When /3i > 2 then the largest contribution is from 

the lower limit of the integration. Thus, when z“ax ^ <C then the leading term in the expansion of the 
approximate expression for the PSD in the power series of w is 


S{f) 


jl + cCfll-l) ) 
1 


0 < /3i < 2 , 
/3i >2. 


This expressions for PSD can also be written as 


Here 


S{f) 



1 — a < j3 < 1 + a, 
/5 1 T a . 


P 


1 + a{Pi 


1 ) = 1 + 


a(i^ — 3) 
2(r7-l) 


(56) 


(57) 


(58) 


is the power-law exponent of the PSD. Equation (1581) generalizes the expression for the power-law exponent obtained 
for nonlinear SDEs [5^. When i/ = 3 then from Eq. (1551) follows that we obtain 1// spectrum. 


A. Power spectral density from scaling properties 


Power-law exponent (l58l) in the PSD can be obtained from the scaling properties of Eq. (|49l), similarly as it has 
been done for the nonlinear SDEs [s^. Changing the variable x to the scaled variable Xg = ax in Eq. (H5|) yields 


|p(x„t) = 


,2(r,-l) OA 


{ (5 - ^ *)] + } . (59) 


The Riemann-Liouville fractional derivative has the following scaling property: qD^ “/(ct) = Thus, 

changing the time t to the scaled time tg = a^t we get 




{(U9) U [.“"-PCx.t.)] + IX . 


(60) 


The change of the variable x to the scaled variable ax or the change of the time t to the scaled time a^t produce the 
same fractional Eokker-Planck equation if 




fJ- = 


a 


(61) 
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It follows, that the transition probability P(a;, tjcco, 0) has the following scaling property: 

aP(ax, t\axo, 0) = P{x, a^t\xo, 0). (62) 

As has been shown in Ref. the power-law steady state PDF Po{x) ^ x~’^ and the scaling property of the transition 
probability (ED) lead to the power-law form PSD S{f) ^ f ^ in a wide range of frequencies. The power-law exponent 
/3 is given by 


P = l + („-3)/^, (63) 

Using Eq. EO we obtain the same expression for (3 as in Eq. (j58p . 

The presence of restrictions at a: = Xmin and x = a^max makes the scaling (|62ll not exact. This limits the power-law 
part of the PSD to a finite range of frequencies /min / <C /max- Similarly as in Ref. [s^, we estimate the limiting 
frequencies as 


a-<^x^l2 27r/ < tJo^mix ?? > 1, (64) 

<C 27r/ < , r] <1. 

This equation shows that the frequency range grows with decrease of a. By increasing the ratio a^max/a^min one can 
get an arbitrarily wide range of the frequencies where the PSD has 1//^ behavior. 


V. NUMERICAL APPROACH 
A. Numerical approximation of sample paths 

Since analytical solution of time-fractional Fokker-Planck equation can be obtained only in separate cases, there 
is a need of numerical solution. Numerical solution of time-fractional Fokker-Planck equation is complicated [H^- 
It is easier to numerically solve Langevin equations E]), (jd]) instead. The desired properties of the solution of the 
Fokker-Planck equation then can be calculated by averaging over many sample paths obtained by solving the Langevin 
equations. The numerical method of solution of the Langevin equations with constant drift coefficient is outlined in 
[^1^. We can use the same method also when the drift coefficient is position-dependent. 

Choosing the time step At of the operational time r the inverse subordinator S{t) is approximated as 

SArit) = [min{n G N : T(nAT) > t} — I]At . (65) 


Such approximation satisfies [1^ 


sup [SArit) - S(t)] ^ At . 

0<t<T 


( 66 ) 


The values T(n At) are generated by summing up the independent and stationary increments of the Levy process: 

r(nAT) =r([n- 1]At)-1- At1/“^„. (67) 


Here are independent totally skewed positive a-stable random variables with the distribution specified by the 
Laplace transform (e“^^) = . Such variables can be generated using the formula 




sin [g {U + f)] 
cos(t/)i 


l^ cos [U-a(l7+|)] ^ 



( 68 ) 


Here U is uniformly distributed on (—§, f) and W has an exponential distribution with mean I. Note, that in Ref. [s^ 
incorrect formula for generating totally skewed positive a-stable random variables has been used. The definition of 
the Levy a-stable distribution using the Laplace transform El differs from the more common definition using the 
Fourier transform. This has been corrected in Ref. [^ . 

The SDE El ia operational time t can be numerically solved using the Euler-Maruyama scheme with the 
time step At. Eor each value of the stochastic variable Xk we assign the physical time tk generated by the process 
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T(r) using Eq. (l67l) . Thus the numerical method of solution of Langevin equations ([2]), (|4]) is given by the following 
equations: 

Xk+i =Xk + a(a;fe)Ar + b{xk)'/^ek , 
tfe+i = tk + . 

Here £k are i.i.d. random variables having standard normal distribution. 

For numerical solution of nonlinear equations, such as those resulting in Eq. 69]), the fixed time step At can be 
inefficient. For example, in Eq. 63) with f] > 1 large values of stochastic variable x lead to large coefficients and 
thus require a very small time step. A more efficient way of solution is to use a variable time step that adapts to 
the coefficients in the equation. Similar method has been used in Refs. [13, for solving nonlinear SDEs. Such a 
variable time step is equivalent to changing of the operational time r to the position-dependent operational time t' . 
If we choose the intensity of random time in Eq. (1^51) as g{x) = b(x)~'£ then, according to Eq. (IMll instead of initial 
Langevin equations Q, 0 we get the new Langevin equations 

^ 

dt{T') = b{x{T'))-^dL‘^{T'). (72) 

Discretizing the operational time r' with the time step At' and using the Euler-Maruyama approximation for Eq. 6ID 
instead of Eqs. 63 we have 

a;fc+i = ajfc-I-^^^^A t'- b v^Ar'efe , (73) 

_1 

( 74 ) 

Comparison with Eqs. (113, (63 shows that Eqs. (63 j (TTH) can be obtained by replacing the time step At in Eqs. (63, 

(63 by 


(69) 

(70) 


At —>■ 


At' 

b{xkY ' 


(75) 


As an example, we solve the Langevin equations 

dx = dr + x'^dW {t) , (76) 

dt = dL°‘{T) (77) 

resulting in the time-fractional Fokker-Planck equation (63- For restriction of the diffusion region we use the reflective 
boundaries at x = Xmin and Xmax- More effective numerical solution scheme is obtained changing the operational time 
T to the time t' defined by the equation 


dt(T') = x(t') ^^dL“(T'). (78) 

This change is equivalent to the introduction of the variable time step At^ = At'x^ . Discretizing the operational 
time t' with the step At' from Eqs. (63“(63 got tbe following numerical approximation: 

Xfe+i = Xk + (v - Xk^T + XkVAr'ek , (79) 

tk+l =tk+ ^ 2 (> 7 -l) ^ ■ (^ 0 ) 

Sample path obtained using Eqs. (63; (1501) with the parameters 77 = 2 and = 3 is shown in Fig.[T| The change of 
the operational time t' with the physical time t is shown in Fig. [TJa) and the dependence of the stochastic variable 
X on the physical time t is shown in Fig. [T](b). Due to nonlinear coefficients in Eq. (l76l) the sample path in Fig. [ijb) 
exhibits peaks or bursts, corresponding to the large deviations of the variable x. The intervals with x being constant 
indicate the heavy-tailed trapping times. Comparing Fig. 0a) with Fig. 0b) we see that the operational time t' 
increases faster when x acquires larger values, in accordance to Eq. 63- 
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t t 

FIG. 1. Sample path obtained from Langevin equations ([761). ([77ll using numerical solution scheme given by Eqs. ([791). ([soil, 
(a) Dependence of the operational time r', defined by Eq. dzll), on the physical time t. (b). Dependence of the stochastic 
variable x on the physical time t. The parameters are a = 0.7, rj = 2, u = 3. Reflective boundaries are placed at = 1 and 

a:max = 1000. 


B. Power spectral density 

Since the equations exhibit a slow (power-law instead of a usual exponential) relaxation , calculation of the PSD 
using sample paths is very slow. More efficient way is to find the eigenvalues and eigenfunctions of the Fokker-Planck 
operator ([5]) and calculate the PSD using the rapidly converging series in Eq. (ITTl) . This is the approach for calculating 
the PSD used in Ref. [s^ for the case of constant diffusion coefficient. 

As an example let us calculate the PSD of the diffusion described by the time-fractional Fokker-Planck equation 
(l4^ with 77 ^ 1 and the reflective boundaries at a:min = 1 and Xmax = The equation (1^ for the eigenfunctions of 
the Fokker-Planck operator that enters Eq. (HU) is 

■ ■ i) 

The reflective boundaries lead to the conditions 5'a(1) = 0 and S\{^) = 0, where 

Sx{x) = - 0 x^’^-^P^{x) - (82) 

is the probability current related to the eigenfunction Px {x). The steady state solution of Eq. is 


■ (83) 

It is more convenient to transform Eq. (ISTl) into the Schrodinger equation [i^ . To do this we first make the diffusion 
coefficient constant by changing the variable x to 




(84) 


Eq. m then becomes 


v' d I 15^ 


with the reflective boundaries at Zmin where 


-^min 


I 1-2) ’ 


r? > 1, 
r? < 1, 


'^max 


^<1- 


(85) 


( 86 ) 
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FIG. 2. (Color online) Dependence of numerically obtained first moments of the variable x on the eigenvalues A for the lowest 
eigenvalues (red dots). Eigenvalues and eigenfunctions are obtained numerically solving Eq. (I88|l . The dashed green line shows 
the slope predicted by Eq. (ISTI) . The parameters used are r? = |, = 3, Xui\u = 1 and Xmax = 1000. 


Here 

T]-l 

Eq. (|^ can be transformed into the Schrodinger equation [i^ 

- + V{z)i^x{z) = \^x{z) 

with the potential 


H(z) 



'(2 + ^ 0 - 


(87) 


( 88 ) 


(89) 


Here V’a(2) = Px{z)/x/P q{z). The condition of zero probability current at the reflective boundaries 0 
2 : = Zmax become 


/_d ^1\ 

\ dz 2 z / 




= 0 . 

^ —2^min )-2'max 


The solution of Eq. (l88l) corresponding to the eigenvalue A = 0 is 


Zinin and 


(90) 




v' - 1 


l-v' _ 1-u 

2v^in ^max 




(91) 


Eq. (1881) can be solved using standard finite-difference or finite-element methods. Having the eigenfunction '(/>a(z) the 
first moment of the stochastic variable x can be calculated using the equation 


^A 


ipQ{z)\ri — 1| i-'i zi-'i il)x{z) dz . 


(92) 


Let us take the following values of the parameters in Eq. 69]); r; = |, V = ?>. The dependence of the numerically 
calculated first moment Xx on the eigenvalue A for lowest eigenvalues is shown in Fig. [2j We see a good agreement 
with the analytical prediction (1511) of power-law dependence on A. For larger eigenvalues A than those shown in Fig. [2] 
the power-law dependence does not hold and Xx decrease faster. 

The PSD calculated using Eq. (1471) is presented in Fig.[3l Eigenvalues A and the first moments Xx shown in Fig. © 
have been used. We see a good agreement with the predicted power-law dependency of the PSD on the frequency for 
frequencies / > /min ~ 1- The power-law exponent coincides with Eq. (I58p . For smaller frequencies / < 1 the PSD 
exhibits the power-law behavior (1481) with the exponent 1 — a. 
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S(f) 



f 


FIG. 3. (Color online) Power spectral density for the diffusion process defined by Eq. (I49II with the parameter a = 0.8. The 
solid red line shows the result of numerical calculation using Eq. ([471). The dashed green line shows the slope 1//, whereas the 
dotted blue line shows the slope Other parameters are the same as in Fig. [2] 


VI. CONCLUSIONS 


In summary, we proposed Eq. m describing the subdiffusion of particles in an inhomogeneous medium that 
generalizes the previously obtained time-fractional Fokker-Planck equation with the position-independent diffusion 
coefficient. Fokker-Planck equation with the position-independent diffusion coefficient has been used to model various 
phenomena such as ion channel gating and the translocation dynamics of a polymer chain threaded through a 
nanopore [b^. Properties of such equations has been studied extensively. In this paper we analyzed a more general 
case when both drift and diffusion coefficients are position-dependent. We hope that the present model can serve as 
a basis to study trapping induced subdiffusion in complex inhomogeneous media. 

We derived the analytical expression of power spectral density of signals described by the one-dimensional time 
fractional Fokker-Planck equation in a more general case when diffusion coefficient depends on the position. The 
general expression for the PSD (I47p we applied to a particular case (1491) when the drift and diffusion coefficients 
have power-law dependence on the position. The resulting PSD has a power-law form S{f) ^ f~^ in a wide range 
of frequencies, with the power-law exponent /3 given by Eq. (I58p . This approximate results is confirmed by the 
numerical simulation (see Fig. ISl). Thus, according to Eq. (1551) . time-fractional Fokker-Planck equation with power- 
law coefficients yields the PSD with the power-law exponent equal to or larger than 1 in a wide range of intermediate 
frequencies. In contrast, the PSD for small frequencies has a power-law dependency on the frequency in the form of 
even when the diffusion coefficient depends on the position. 

Since an analytical solution of time-fractional Fokker-Planck equation can be obtained only in separate cases, there 
is a need of numerical solution. For the numerical solution of the nonlinear equations, such as those resulting in 
Eq. (14^ . we propose to use a variable time step that adapts to the coefficients in the equation. Such a variable time 
step is equivalent to changing of the operational time r to the position-dependent operational time t' . 
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